Finite-temperature dynamics with the density-matrix renormalization group method 
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We present a new numerical method for the evaluation of dynamical response functions at finite temperatures 
in one-dimensional strongly correlated systems. The approach is based on the density-matrix renormalization 
group method, combined with the finite-temperature Lanczos diagonalization. The feasibility of the method is 
tested on the example of dynamical spin correlations in the anisotropic Heisenberg chain, in particular it yields 
nontrivial results for the critical behavior in the isotropic case. 
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I. INTRODUCTION 

Strongly correlated systems present one of the major theo- 
retical challenges in last decades and are stimulating the inten- 
sive search for adequate numerical methods to evaluate their 
properties. Within the low-dimensional systems, in particu- 
lar one-dimensional (ID) systems the breakthrough has been 
achieved with the introduction of the density matrix renormal- 
ization group (DMRG) method' allowing accurate calcula- 
tion of the ground-state wavefunction and its static properties 
on large systems far beyond those available with the exact- 
diagonalization methods. Among various DMRG extensions^ 
we concentrate here on the goal to study the dynamical re- 
sponse of such systems at finite temperatures T > 0. It should 
be observed that in spite of the satisfactory description and 
an understanding of static properties of generic ID systems 
at r > the corresponding dynamics, in particular the low- 
frequency one as manifested in the transport quantities, NMR 
relaxation, is far less understood and approachable via numer- 
ical methods. 

For dynamical response within the ground-state the tar- 
geting within the DMRG has been extended to contain also 
excited states.-' "* Transfer-matrix DMRG^ *"'^ is very efficient 
to evaluate thermodynamic properties of models with short- 
range interactions, as well as some dynamical correlations 
of very limited range. Time dependent DMRG**'^ developed 
recently enables studies of short-time evolution of general 
many-body systems, hence also of T > behavior, but is 
rather limited in reaching the low-cj response. Recently, a 
DMRG method extended with the polynomial expansion has 
been proposed to treat low-T dynamics.'" On the other hand, 
methods emerging from the exact diagonalization approach 
as the T > Lanczos method (FTLM)" and the low-T 
version'^ have high-w resolution and provide the information 
on the nontrivial dynamics of correlated models, but are still 
restricted to small systems reachable with exact diagonaliza- 
tion. 

The paper is organized as follows. In the next section (HJl 
we first present our new method, with which we calculate 
some static and dynamical properties of the model described 
in sectionHni In section |lll] we also show our results, first the 
test of our method on the XY model and then our main results 
for the isotropic Heisenberg model. In the last section HVl we 



present our conclusions. 



n. METHOD 

In this paper we propose a new method for the calcula- 
tion of the r > dynamics which is a combination of the 
FTLM and the DMRG, namely the finite-temperature dynam- 
ical DMRG (FTD-DMRG) method. It is constructed to cal- 
culate dynamical response functions in ID systems at T > 0, 
with the emphasis on the low-cj regime. As a test we con- 
sider highly nontrivial spin correlations within the anisotropic 
Heisenberg model on a chain. 

In the standard T = DMRG the ground-state is used to 
construct the basis. In our case we use the full T > density 
matrix, which can in general be expressed with eigenstates \n) 
and corresponding eigenvalues i?„. 
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where (3 — 1/T and Z is the (grand)canonical sum. We pro- 
ceed by extending the density matrix, Eq.([Tl), with the sam- 
pling over the random vectors |r) = J2n f^rnln) where f3rn 
denote random amplitudes. 
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It is easy to show that Eq.(l2]l reduces to Eq.([TJ expressed in 
diagonal basis |?^) (n| since offdiagonal terms vanish assuming 
normalized and random \r)^ 

In Eq. (|2]i we evaluate the operator e"^^!"^ on |r) by start- 
ing the Lanczos procedure from |r). After diagonalization of 
the Lanczos tridiagonal iJ, we obtain the first series of Lanc- 
zos eigenvectors with corresponding eigenenergies e[, 
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It is evident that for M approaching Nst Eq. (O reproduces 
fully Eq. ([T]i, while for M <C A^^f as used in practice rep- 
resents an efficient way of evaluation of density matrix. The 
sum Z may be evaluated in the same manner as within the 
FTLMii 
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In the original T = DMRG procedure one targets the 
ground-state.ii^ Instead, at T > we target states \ipr) and 
construct the density matrix according to Eq. (O. 

Since our aim is to calculate dynamical response functions, 
expressed as autocorrelation functions, we also requke a good 
representation of the operator density matrix. 
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It replaces the operator on the ground-state in original T — 
DMRG^'^ and is evaluated by extending Eq. (|3]l. 
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In the proposed targeting we sum up above contributions with 
weighting factors. 
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with the restriction pi + P2 = 1. From ptot we prepare the 
reduced density matrix by integrating out the environment, 
which is then used to construct the basis within the infinite 
and finite algorithms of the DMRG.^ Our way of targeting 
is in fact very similar to the one in Ref. [l^ with an addi- 
tional random sampling suppressing the non-diagonal terms 
of p. In such a way we prepare the basis for any T > 0, 
whereby limitations are emerging from the truncation of the 
basis being more under control for low T. It should also be 
mentioned that for dynamical response at particular w there is 
an improvement to target also excited states corresponding to 
so called correction vectors. Still, the latter does not affect 
quality of the most interesting and challenging regime uj ^ 
as well as it increases the computational demand, hence we do 
not employ it here. 

Physical quantities are calculated in the measurement part 
of the FTD-DMRG procedure in the same manner as within 
the FTLM. ' ' A dynamical autocorrelation function 
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is evaluated with the use of two Lanczos series of eigenstates 



and eigenenergies. 
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The second Lanczos series of eigenstates IV'/'^) and eigenen- 
ergies e^"^ is obtained from second Lanczos procedure starting 
from the initial vector A\r). 



III. MODEL AND RESULTS 

As a nontrivial test of the method we analyse the dynamics 
of the ID anisotropic Heisenberg model. 
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where Sf,S^ are local spin S = 1/2 operators, L is the 
chain length, J is the exchange coupling (in the following we 
use J = 1) and A the anisotropy parameter In our calcula- 
tions we focus on systems in the absence of the magnetic field, 
hence on the subspace S^^f = 0. As the quantity of interest 
we choose the dynamical spin structure factor S{q, io) and the 
corresponding susceptibility xili 
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dte^"*(5^(g,t)5^((Z,0)), 



X"{q,cj) = n{l-e-''nS{q,u;). 
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As usual within the DMRG technique more accurate re- 
sults are obtained with open boundary conditions,^ hence one 
defines S^{q) = y/2/{L+l)J2sin{qi)Si whereby q = 
ttJ /{L + 1) with j — 1, . . . L. In our calculations we con- 
centrate on most challenging q = Q = tt, i.e. j — L. 

The relaxation function $(g,w) = x"{1t^)I^ should be 
an even function of cj. This represents another nontrivial test 
for the FTD-DMRG method. In addition to considering com- 
plete spectra ~x"{q,Lo) better defined criteria are frequency 
moments. 
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Due to symmetry only even M" {q) are finite while the static 
susceptibility corresponds to x°('Z) = M^{q). 

In the following we employ the FTD-DMRG method to 
evaluate for A = 0, 1 and various T. In the ac- 

tual implementation we use the infinite and finite-size DMRG 
basis preparation and the calculation of S{Q,iS) via Eq.® 
{A ~ Sz{Q)) performed on the system divided into two sub- 
blocks of size (i — 2)/2 and two coupling sites in betweeni^ 
In the preparatory sweeping typically 1 or 2 sweeps are suffi- 
cient for the convergence of the basis. Important parameters 
for the final quality of results are the (subblock) DMRG trun- 
cation number m and the number of Lanczos steps M . We 
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are typically restricted to m ~ 1000 and M ~ 100. We have 
two kinds of sampling. One in the determination of the den- 
sity matrix Eq. (|2]i, R — Ri for the basis preparation, and the 
other in the evaluation of the final Eq. (|9]l, R = R-z- While 
only modest Ri ^ 50 is adequate, i?2 ^ 1 is needed in par- 
ticular at low T— to get accurate matrix elements. At higher 
T R2 can be reduced effectively to _R2 ^ Furtheron we 
mainly consider T < 0.5 with R2 ^ 100. When evaluating 
the feasibility of various methods we should keep in mind that 
the full exact diagonalization evaluation of S{q, cj) at T > 
for the model at hand can be performed up to L = 14 — 16, 
with the FTLM technique up to i = 24, while in the follow- 
ing we present the FTD-DMRG results up to i = 40. 



A. XY model 

The A = case maps onto noninteracting spinless 
fermions via the Jordan- Wigner transformation and S{Q,uj) 
can be expressed for any T > in a Lindhard form. For finite 
L the only caveat is that the FTD-DMRG is performed within 
a canonical systems with fixed S^^^ = 0, i.e., with fixed num- 
ber of fermions ~ L/2 while the usual (easier) evaluation 
is within the grandcanonical ensemble. In Fig. [T]we present 
the FTD-DMRG result for (unsymmetrized) relaxation func- 
tion $((5, Lu) at low T = 0.25. Results are for L = 36 where 
the basis is heavily reduced, i.e., only 5 x 10~^ of all states 
are retained within the final evaluation. For comparison we 
show the exact (grandcanonical) fermionic result for the same 
system with open boundary condition and for all presented 
spectra we use the damping 77 = 0.05. Oscillations are a clear 
sign of finite-size system and slowly disappear with increasing 
T and L. The finite-size effect can be avoided by smoothing 
with a Gaussian filter with the width adapted to the frequency 
oc 1/L. From Fig. [T]it is evident that at low \uj\ < 1 the 
agreement between the FTD-DMRG and the exact result is 
very satisfactory. At high uj ^ 2 the FTD-DMRG does not 
fully reproduce the sharp spectral edge which could be im- 
proved by the introduction of the correction-vector targeting 
for w 7^ within the method. 

Fig. |2] shows the corresponding results for the frequency 
moments Mn (Q) displayed vs. 1/L obtained with the full ba- 
sis for i < 22 and with the FTD-DMRG method for L < 40. 
For comparison also corresponding exact results are shown 
within the canonical calculation at — L/2. It is evident 
that T = 0.25 is already high enough so that moments are 
essentially size independent. Also up to i = 40 FTD-DMRG 
results are well stable, at least for lowest Mq, M2, while for 
M4 some deviations originate from high-w regime and are 
also visible in Fig. [T] At the same time. Mi, M3 « is well 
reproduced as required by the symmetry of ^{Q, oj). 



B. Isotropic Heisenberg model 

The isotropic A = 1 case (at S^^f = 0) representing 
marginally gapless system is by far more challenging. For 
T > there are no exact results for dynamical quanti- 
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Figure 1: (color online) Relaxation function ^{Q,uj) within the XY 
model for T — 0.25 and a system of L = 36 sites. For comparison 
the exact grandcanonical result for spinless fermions is shown and 
the corresponding smoothed curve relevant for L ^ 00. 
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Figure 2: (color online) Lowest frequency moments Mn{Q) vs. 1/L 
for T — 0.25. For comparison exact moments are shown up to L = 
22. 

ties. The bosonization approach provides a form for S{q, u) 
within the low oj — T regime .-^^^ Relative to the A = 
case the divergence for A = 1 is stronger and nontrivial. 
The isotropic model has been an obvious target for numeri- 
cal methods. Static quantities, as the structure factor S{q) and 
X°(g) have been evaluated using the quantum Monte Carlo 
(QMC) method and the high-T expansion,'"' '^ recently also 
with the time-dependent DMRG,'^ but only for q ^ Q so 
far An obvious deficiency is in results for dynamic quanti- 
ties at w ~ since the QMC approach (due to the Maximum 
Entropy procedure) seems to have considerable uncertainty 
in this regime:'^ On the other hand, the latter regime is fre- 
quently just the most interesting, e.g., in connection with the 
NMR relaxation rate 1 /Ti cx J2q ^g^iq, 0), with trans- 
port quantities etc. 

In Fig. [3] we present results for ^{Q,uj) obtained for 
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L = 40 sites and different T. Since spectra are peaked at 
w = (in contrast to Fig. [TJ finite-size oscillations are more 
pronounced. Hence, also smoothed curves (Gaussian width 
a = 4cos(7rL/2(L + 1))/V2) are presented as relevant for 
L oo. We note that such spectra are nearly i-independent 
(L = 16 — 40) for lli > 0.5 whereas for u ^ stiU scale as 
a + b/L. 




Figure 3: (color online) Symmetrized ${Q,u!) for the isotropic 
Heisenberg model shown for L = 40 and T = 0.15,0.25,0.35. 
Presented are also finite-size smoothed spectra (dotted line). 



On the other hand, static x"((5) can be extracted directly 
without invoking any smoothing and FTD-DMRG results 
combined with the FTLM results for L — 12 — 20 are shown 
in Fig. |4] scaled vs. 1/L. Deviations from the linear scaling 
mostly emerge from the random sampling in the basis prepa- 
ration and the dynamical quantity evaluation, and for the latter 
are indicated with error bars. Final scaled FTD-DMRG results 
for x^iQ) vs. T are shown in Fig.|5l together with the result 
of the QMC analysis— of the analytical expression 
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Figure 4: (color online) X°(Q) for ths isotropic model vs. 1/L for 
different T as calculated via the FTD-DMRG method for L = 22 - 
40 and via FTLM for L = 12 - 20. 
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Figure 5: (color online) Scaled values for x'^iQ) <*nd S{Q,uj — 0) 
vs. T for the isotropic model. The dashed curve represents 
using the analytical form as extracted from the QMC results Ref . [Ij- 
Dotted S{Q, uj = Q) curve is the fit as deduced from the analytical 
approximation. 



Our FTD-DMRG result is quite consistent with QMC results 
at higher T > 0.3. Still it is indicative that we get higher 
values (beyond error bars) for T < 0.3. 

Finally, we present in the same Fig. |5] also scaled values of 
S{Q, UJ — 0) vs. T. Bosonization theory givesii 

5(Q,0) = ^[ln(A/T)]i/2 (14) 

also fitted to our results with A = 24.27 taken from Ref. [H 
and adjusted A ^ 0.205. The agreement with the analytical fit 
is very good although there seems to be substantial difference 
in the prefactor A.^* On the other hand, it should be reminded 
that for this quantity there are no reliable larger-system al- 
ternative results since the QMC analysi a'^'^^ appears to have 
some difficulties in the regime oj ^ 0. 



IV. CONCLUSIONS 

In conclusion, we have introduced the FTD-DMRG 
method, which is the extension of the density matrix-based 
optimization of target states and the FTLM method for the 
evaluation of dynamical quantities at T > 0. It is so far well 
founded and tested for relatively low T and not too large sys- 
tems, e.g., L < 40, while the feasibility or possible break- 
down at larger T should still be understood. Presented results 
are obtained for systems with Z < 200 (normalized so that 
Z{T = 0) = 1) although the method is not in principle lim- 
ited to low T since it is not essential that all relevant many- 
body states are well represented, in analogy to the FTLM.— 
The emphasis so far is on the most challenging oj ^ dynam- 
ical response while higher uj could be improved by extending 
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the density matrix by optimizing the correction vector at par- 
ticular wiii!^ As the test we use the A = case which is 
nontrivial for the FTD-DMRG method while exact results are 
available via the spinless-fermion representation. On the other 
hand, results for the isotropic A = 1 case where we concen- 
trate on the low lo ~ T regime of dynamical spin correlations 
S'(Q, Lo) show that the presented method goes beyond the ca- 
pabilities of up-to-date numerical methods, e.g., in the case of 
S{Q,uj ~ 0). Clearly, more effort is needed to examine in 
more detail the feasibility of the new method. 
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